library(ca)
results <- list.files(pattern=".*results.txt")
tawc <- vector("list",length(results))
for (r in results){ 
	tawc[[which(results==r)]] <- read.delim(r,header=T,sep="\t")
	tawc[[which(results==r)]]$filename <- r
}
tawc <- do.call("rbind",tawc)

tawc <- tawc[,c("col0.Competence",
"col0.Excitement",
"col0.Rugged",
"col0.Sincerity",
"col0.Sophistication")]
colnames(tawc) <- gsub("col0.","",colnames(tawc))
rownames(tawc) <- c("Cons","Green","Lab","LDems","UKIP")
tawc <- tawc[c("Cons","Lab","LDem","UKIP","Green"),]

ca.tawc <- ca(tawc)

###  Bootstrapping
### Code taken from http://www.mattpeeples.net/caboot.html
rown <- nrow(tawc) # number of rows in data table
coln <- ncol(tawc) # number of columns in data table
nsim <- 1000 # number of bootstrapped simulations to run
cut <- 0.95 # set convex hull confidence interval for final plot
lab <- colnames(tawc) # define column (type) labels

# create nsim bootstrapped replicates from original data
data.rowsum <- apply(tawc,1,sum)
data.sim <- rmultinom(nsim,data.rowsum[1],prob=tawc[1,]) # create bootstrapped replicate for first row with probability defined by actual proportions of types
for (i in 2:rown) {
data.sim <- cbind(data.sim,rmultinom(nsim,data.rowsum[i],prob=tawc[i,]))} # create bootstrapped replicates for remaining rows

# reorder bootstrapped replicates by number of rows
data.sim <- t(data.sim)
data.sim2 <- matrix(rep(0,nsim*rown*coln),nrow=nsim*rown)
for (k in 1:nsim) {
	for (i in 1:rown) {
		data.sim2[(k-1)*rown+i,] <- data.sim[k+(i-1)*nsim,]}}

# define simulated column principal coordinates using row transition formula
data.rowsc <- ca.tawc$rowcoord[,1:2]
data.colsim <- t(t(data.rowsc) %*% data.sim2[1:rown,]) / apply(data.sim2[1:rown,],2,sum)
for (k in 2:nsim) {
data.colsim <- rbind(data.colsim, t(t(data.rowsc) %*% data.sim2[((k-1)*rown+1):(k*rown),]) / apply(data.sim2[((k-1)*rown+1):(k*rown),],2,sum))}
data.colsim2 <- matrix(rep(0,nsim*coln*2),nrow=nsim*coln)

# define simulated row principal coordinates using column transition formula
data.colsc <- ca.tawc$colcoord[,1:2]
data.rowsim <- t(t(data.colsc) %*% t(data.sim2[1:rown,])) / apply(data.sim2[1:rown,],1,sum)
for (k in 2:nsim) {
x <- t(t(data.colsc) %*% t(data.sim2[((k-1)*rown+1):(k*rown),]))
data.rowsim <- rbind(data.rowsim, x / apply(data.sim2[((k-1)*rown+1):(k*rown),],1,sum))}
data.rowsim2 <- matrix(rep(0,nsim*rown*2),nrow=nsim*rown)

# populate row and column simulated point matrices with projected coordinates
for (j in 1:coln) {
	for (k in 1:nsim) {
		data.colsim2[(j-1)*nsim+k,] <- data.colsim[j+(k-1)*coln,]}}
for (j in 1:rown) {
	for (k in 1:nsim) {
		data.rowsim2[(j-1)*nsim+k,] <- data.rowsim[j+(k-1)*rown,]}}


# set up matrices for recording convex hull sizes
x2 <- matrix(0,nrow(tawc),1)
x3 <- matrix(0,nrow(tawc),1)
x2_a <- matrix(0,nrow(x2),1)
x3_a <- matrix(0,nrow(x3),1)
row.names(x2) <- row.names(tawc)
row.names(x3) <- row.names(tawc)
row.names(x2_a) <- row.names(x2)
row.names(x3_a) <- row.names(x3)

# Define convex hulls around row points and record hull sizes
poly.row <- list()
for (j in 1:nrow(tawc)) {
	points <- data.rowsim2[(nsim*(j-1)+1):(nsim*j),]
	hpts <- chull(points)
	repeat {
		hpts <- chull(points)
		npts <- nrow(points[-hpts,])
		if(npts/nsim<cut) break
		points <- points[-hpts,]}

	hpts <- c(hpts,hpts[1])
	v <- points[hpts,]
	poly.row[[j]] <- v
	v <- rbind(unique(v),v[1,])
	v2 <- as.matrix(dist(v))
	x <- matrix(0,nrow(v2)-1,1)
	for (i in 1:nrow(x)) {
	x[i,1] <- v2[i+1,i]}
	x2[j,1] <- sum(x)
	x3[j,1] <- max(v2)}

# Define convex hulls around column points
poly.col <- list()
for (j in 1:coln) {
	points <- data.colsim2[(nsim*(j-1)+1):(nsim*j),]
	hpts <- chull(points)
	repeat {
		hpts <- chull(points)
		npts <- nrow(points[-hpts,])
		if(npts/nsim<cut) break
		points <- points[-hpts,]}

	hpts <- c(hpts,hpts[1])
	poly.col[[j]] <- points[hpts,]}


### Plotting

png(file="party_pos.png",width=1800, height= 1000, res = 300)
par (mar=c(3,3,1,1), mgp=c(2,.7,0), tck=-.01)


plot(data.rowsc[,1],data.rowsc[,2],type='n',
	xlab="Primary dimension",ylab="Secondary dimension",
	main="",
	xlim = c(-.3,.3),
	ylim = c(-.3,.3))
data.col <- t(t(data.rowsc) %*% as.matrix(tawc)) / apply(tawc,2,sum) # set column label positions

party.colours <- c("#0087DC",
	"#D50000",
	"#FDBB30",
	"#70147A",
	"#99CC33")
party.ltys <- 1:5
party.names <- rownames(tawc)
fade.colour <- function(col,alpha){
	components <- col2rgb(col)
	return(rgb(components[1,]/255,components[2,]/255,components[3,]/255,alpha=alpha))
	}

### Plot the convex hulls
for (i in 1:length(poly.row)) { 
	polygon(poly.row[[i]],
		## col=fade.colour(party.colours[i],0.5),
		lty = party.ltys[i],
		lwd = 2,
		border=party.colours[i])
} 
legend("bottomright",
	legend = party.names,
	lty = party.ltys,
	lwd = 2,
	col = party.colours,
	bty = "n")
	
## Row label positions
centroids <- lapply(poly.row,colMeans)
names(centroids) <- rownames(tawc)
for (i in 1:length(centroids)) {
	text(x= centroids[[i]][1],
		y = centroids[[i]][2],
		label = names(centroids)[i],
		cex = 0.5)
}

dev.off()

### Now traits
png(file="trait_pos.png",width=1800, height= 1000, res = 300)
par (mar=c(3,3,1,1), mgp=c(2,.7,0), tck=-.01)

plot(data.rowsc[,1],data.rowsc[,2],type='n',
	xlab="Primary dimension",ylab="Secondary dimension",
	main="",
	xlim = c(-.3,.3),
	ylim = c(-.3,.3))

for (i in 1:length(poly.col)) { 
polygon(poly.col[[i]],lty = 1)
## Column label positions
text(data.col[i,1],data.col[i,2],lab[i], 
	cex = 0.5)
}
for (i in 1:length(centroids)) {
	text(x= centroids[[i]][1],
		y = centroids[[i]][2],
		label = paste0("(",names(centroids)[i],")"),
		font = 3,
		cex = 0.5)
}

dev.off()
